function [PI,M,indprob] = my_solab(A,B,NK,stake)

% Based on codes by Paul Klein (solab.m) and Paul Soederlind (ComItAlg.m)
% Inputs: Two square matrices A and B:
%
% A*x(t+1) = B*x(t)
% 
% where x(t) is arranged so that the state variables come first, and
% nk is the number of state variables.
%
% Outputs: the decision rule PI and the law of motion M. If we write
% x(t) = [k(t);u(t)] where k(t) contains precisely the state variables:
% u(t)   = PI*k(t) 
% k(t+1) = M*k(t)
% indprob = 1 if eigenvalues are not implying a saddlepath solution in the
% sense of Blanchard-Kahn

indprob = 0;


%Complex Schur Decomposition
[s,t,q,z] = qz(A,B);   

%Pick non-explosive (stable) eigenvalues
slt = (abs(diag(t))<stake*abs(diag(s)));  
nk=sum(slt);

%Reorder the system with stable eigs in upper-left
[s,t,q,z] = ordqz(s,t,q,z,slt);   

%Split up the results appropriately
z21 = z(nk+1:end,1:nk);
z11 = z(1:nk,1:nk);

s11 = s(1:nk,1:nk);
t11 = t(1:nk,1:nk);

%Identify cases with no/multiple solutions
if nk>NK
    warning('The Equilibrium is Locally Indeterminate')
    exitflag=2;
elseif nk<NK
    warning('No Local Equilibrium Exists')
    exitflag = 0;
end

if rank(z11)<nk;
    warning('Invertibility condition violated')
    exitflag = 3;
end

%Compute the Solution
z11i = z11\eye(nk);
PI = real(z21*z11i);  
M = real(z11*(s11\t11)*z11i);
% 
% [s,t,q,z] = qz(A,B);                      % upper triangular factorization of the matrix pencil b-za
% logconA = abs(diag(t)) <= stake*(abs(diag(s))); % selecting stable eigenvalues
% [s,t,q,z] = ordqz(s,t,q,z,logconA);       % reorder: stable eigenvalues first
% logcon = abs(diag(t)) <= stake*(abs(diag(s)));  % 1 for stable eigenvalue 
% 
% if sum(logcon) < nk;
%   warning('Too few stable roots: no stable solution');
%   M = NaN; PI = NaN;
%   indprob = 1;
%   return;
% elseif sum(logcon) > nk;
%     
%   warning('Too many stable roots: inifite number of stable solutions');
% 
%   M = NaN; PI = NaN;
%   indprob = 1;
%   return;
% end;
% 
% z21 = z(nk+1:end,1:nk);
% z11 = z(1:nk,1:nk);
% 
% if rank(z11)<nk;
% 	error('Invertibility condition violated')
% end
% 
% z11i = z11\eye(nk);
% s11 = s(1:nk,1:nk);
% t11 = t(1:nk,1:nk);
% 
% PI = real(z21*z11i);
% 
% dyn = s11\t11;
% M = real(z11*dyn*z11i);